################
#PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
#
#Observational Analysis ESS
#Table C 5
#
#Verena Fetscher
#July 2022
####################

####################
# Load packages
####################

library(tidyverse)
library(xtable)
library(Zelig)
library(texreg)

####################
# Load data
####################
source("dataframes.R")

##########################
#Model estimation
##########################

##########################
#Fixed effects
##########################
mod1<-zelig(redistribution ~ benefitsAvg+
              factor(year),
            model = 'ls',
            data=mi.out$imp5[mi.out$imp5$income.dist.th>0,], 
            cite = FALSE)
summary(mod1)


mod2<-zelig(redistribution ~ benefitsAvg+
              our+agea+eduyrs+single+gender+uempl+
              factor(year),
            model = 'ls',
            data=mi.out$imp5[mi.out$imp5$income.dist.th>0,],
            cite = FALSE)
summary(mod2)


mod3<-zelig(redistribution ~ benefitsAvg+
              our+agea+eduyrs+single+gender+uempl+
              gini.disp+immigration+social_spending+toptax+red_transfers+
              as.factor(year),
            model = 'ls',
            data=mi.out$imp5[mi.out$imp5$income.dist.th>0,],
            cite = FALSE)
summary(mod3)


texreg(list(mod1,mod2,mod3),stars = c(0.01,0.05,0.1),
       custom.coef.names = c("Intercept", "Benefit concentration","Year 2004",
                             "Year 2006","Year 2008","Year 2010","Year 2012",
                             "Year 2014","Occu. unemployment risk","Age","Education","Single","Female","Unemployment",
                             "Gini (disposable income)", "Immigration",
                             "Social Spending","Top Income Tax",
                             "Redistribution via transfers"))





##########################
#Quantities of interest
##########################

# min to max
mi.out$imp5%>%
  group_by(country)%>%
  summarize(benefits=mean(benefitsAvg))

x_Germany <- setx(mod2, benefitsAvg=0.0861) 
x_UK <- setx(mod2, benefitsAvg=0.777) 

s_out <- sim(mod2, x = x_Germany, x1 = x_UK)
summary(s_out)

# mean to mean+sd
mi.out$imp5%>%
  summarize(mean_benefits=mean(benefitsAvg),
            sd_benefits=sd(benefitsAvg)) -> benefits

benefits

benefits[[2]]
x_mean <- setx(mod2, benefitsAvg=benefits[[1]]) 
x_sd <- setx(mod2, benefitsAvg=(benefits[[1]]+benefits[[2]]))

s_out <- sim(mod2, x = x_mean, x1 = x_sd)
summary(s_out)


